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STOCHASTIC CONTROL SYNTHESIS OF SYSTEMS WITH 
STRUCTURED UNCERTAINTY 


Luis G. Crespo* 


ABSTRACT 

This paper presents a study on the design of robust controllers by using random vari- 
ables to model structured uncertainty for both SISO and MIMO feedback systems. 

Once the parameter uncertainty is prescribed with probability density functions, its ef- 
fects are propagated through the analysis leading to stochastic metrics for the system’s 
output. Control designs that aim for satisfactory performances while guaranteeing ro- 
bust closed loop stability are attained by solving constrained non-linear optimization 
problems in the frequency domain. This approach permits not only to quantify the 
probability of having unstable and unfavorable responses for a particular control design 
but also to search for controls while favoring the values of the parameters with higher 
chance of occurrence. In this manner, robust optimality is achieved while the charac- 
teristic conservatism of conventional robust control methods is eliminated. Examples 
that admit closed form expressions for the probabilistic metrics of the output are used 
to elucidate the nature of the problem at hand and validate the proposed formulations. 

Keywords: Robust control, probabilistic methods, structured uncertainty, robust op- 
timization. 

1 INTRODUCTION 

The main requirement of feedback control is to achieve acceptable levels of performance in 
the presence of uncertainty. Fundamental trade offs and compromises between these two 
aspects motivate the entire body of feedback theory. While performance concerns aspects 
such as reference tracking, disturbance rejection, bounded control effort, etc., uncertainty 
appears as a result of the inevitable discrepancies between the physical problem and its 
deterministic mathematical model. Ignorance on the system’s exact dynamics, on the actual 
operating conditions and the purposeful choice of a simplified representation of the physical 
problem exemplify this aspect. In this context, uncertainty can be classified as structured 
(or parametric ) and unstructured [16]. The first kind corresponds to inaccuracies on the 
parameters of the model while the second one corresponds to unmodeled dynamics. 

Uncertainty can be modeled in many ways depending upon the desired quality of its math- 
ematical description. Differential sensitivity, multi-models, interval analysis, perturbations, 
fuzzy sets and probabilistic methods have been used [2, 7, 21]. The effects of uncertainty 
on the stability associated with the prescribed control solutions have been studied by both 
deterministic [13] and stochastic [8, 17, 18, 19] means. These analysis tools however, have 
not been integrated to the control design process. 

The methods most commonly used for robust control design are p-synthesis [1, 4, 9, 15] 
and i/oo optimization [3, 6, 14, 20]. In these, uncertainty is modeled with norm-bounded 
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complex perturbations of fixed but arbitrary structure about its nominal form. This treat- 
ment is extensively used primarily because it leads to a tractable set of sufficient conditions 
for robust stability. Such approach however, have the following drawbacks: (i) the crudeness 
of the uncertainty description usually leads to redundant and physically impossible plants, 
then to highly conservative designs, (ii) it is not feasible to favor scenarios with higher chance 
of occurrence among all the possible ones, and as a result, robust optimality is precluded, 
(iii) a quantitative description of the robustness of the solution is unattainable and (iv) the 
resulting controllers are so complex that model reduction techniques are usually required. 
While such perturbations account for unstructured uncertainty coarsely, an augmented plant 
model with structured uncertainty can be used to conciliate the uncertainty representation 
with the physics of the problem. While robust optimization has been studied in various dis- 
ciplines using different uncertainty models [5, 10, 11], stochastic control synthesis remains, 
to a large extent, unexplored. 

This paper studies the control design of plants with structured uncertainty for both 
single-input-single-output (SISO) and multiple-input-multiple-output (MIMO) systems us- 
ing a probabilistic approach. The joint probability-density-function (PDF) of the parame- 
ters is prescribed a priori , and then propagated, leading to a probabilistic description of the 
metrics of the controlled response. Control design, involving decoupling, performance and 
stability aspects, is carried out by solving constrained non-linear optimization problems in 
the frequency domain. The content of this paper is organized as follows. Background on the 
performance requirements for feedback SISO and MIMO systems is presented in Section 2. 
Section 3 introduces closed loop stability considerations and the concept of robust stability in 
the probabilistic sense. Several design strategies, based on the minimization of a frequency 
dependent cost function, are proposed in Section 4 and evaluated in Section 5. A discussion 
on the developments and some conclusions, are presented in Sections 6 and 7. 

2 CONTROL REQUIREMENTS 

Consider a linear-time-invariant (LTI) system with transfer function matrix G(s, p), where 
p is a vector of random variables corresponding to the uncertain parameters. Let p be 
prescribed a priori by the joint probability density function /p(p). In this notation, /p(p) 
is the PDF of p at p = p. The support of /p(p), denoted with A, is the set of values the 
random vector p can take, i.e. p G A. Note that the propagation of A leads to a family 
or a set of uncertain plant models in which the physical system is assumed to reside. The 
probability of occurrence of a plant within this set is fully determined by the structure of 
G and the above mentioned PDF. If the components of p are independent, they can be 
prescribed individually. Let K(s, k) be the transfer function matrix of a controller, where 
k is a vector of gains to be determined. Control requirements for SISO and MIMO systems 
are presented next. 

2.1 SISO Systems 

Let L(s 1 p 1 k)=GK be the well-known loop transfer function. For the one degree-of-freedom 
system shown in figure 1, where r, u, z, d and n stand for reference, control, output, distur- 
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bance and noise signals; the following equations apply 


( 1 ) 

( 2 ) 

( 3 ) 


z = Tr + Sd — Tn 
e = z — r — —Sr + Sd — Tn 
u = KS(r — d — n) 

where S=(I + L)~ l is the sensitivity function and T=SL = 1 — S' is the complementary 
sensitivity function. The following analysis can now be made. For e?sO, disturbance rejection 
and tracking are achieved if <S'~0. This implies that |I/|^>1. On the other hand, noise 
rejection implies that S~ 1, which is obtained if Lzz 0. This illustrates a typical conflict in 
the design process. Hence, the control requirements are 

1. Disturbance rejection: \L\ > 1. 

2. Noise attenuation: \L\ < 1 

3. Reference tracking: \L\ > 1 

4. Control energy reduction: \L\ < 1 

5. Robust stability: Re{s*(p, k)} < 0 for all s*(p, k) such that 1 + L(si, p, k) = 0. 

Fortunately, the conflicting objectives are generally in different frequency ranges, and we can 
meet most of them by having large loop gains at low frequencies (usual range of d and r) 
and small gains at high frequencies (usual range of n). 

Due to the probabilistic nature of p, L is a random complex variable parameterized by 
the vector k and the frequency lo. This makes the closed loop poles random variables and 
the Bode and Nyquist diagrams of G and L random processes continuously parameterized. 
In other words, at a fixed k and u, L is a random variable utterly described by a PDF, 
in contrast to the single complex function it is when p is deterministic. In this context, 
stochastic control design intends to shape the random process of the system’s output by 
manipulating k such that the control requirements are satisfied for all the plants generated 
by the propagation of A. 

2.2 MIMO Systems 

The system dynamics are now described by transfer function matrices. Let C(s,c) be the 
transfer function matrix of a compensator. For the system shown in figure 2, where the open 
loop transfer matrix is L(s , p, k) = GCK , Equations (2)-(l) hold once K is replaced by KC. 
Denoting with d(L) and <j(L) the maximum and minimum singular values of L, the control 
requirements are 

1. Disturbance rejection: a(L) > 1 

2. Noise attenuation: d(L) < 1 

3. Reference tracking: a(L ) > 1 

4. Control energy reduction: d(L) < 1 
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5. Robust stability: i?e{sj(p, c, k)} <0 for all s* such that Det{I + L(s*, p, c, k)} = 0. 

As before, the open loop requirements (1) and (3) are valid at low frequencies while require- 
ments (2) and (4) are valid at high frequencies. Then, at frequencies where high gains are 
required the ‘worst-case’ direction is related to g_(L), whereas at frequencies where low gains 
are required the ‘worst-case’ direction is related to a(L). Both ‘worst-case’ directions can be 
unified by the function 

r = p(w)n(L) + (1 - pM)a(L) (4) 

where p(ca) G [0, 1] is a weighting function whose value approaches one at low frequencies 
and zero at high frequencies, e.g. p(co) = 1/(1 + a ;) m for m> 1. In this context, stochastic 
control design intends to shape the random process T (or equivalently the random process of 
the singular values of L ) by manipulating k such that the control requirements are satisfied 
for all the plants generated by the propagation of A. 

Compensator design is another area where a probabilistic treatment of uncertainty is 
valuable. The objective of the compensator is to eliminate the effects of undesired cross- 
couplings among inputs and outputs such that the MIMO system can be treated as a set of 
independent SISO systems. This can be achieved if CG is diagonal, e.g. if G is square and 
and D is the desired uncoupled behavior, C ~ DG~ l . Since the exact form of G is always 
unknown, compensators based on a nominal deterministic plant will not perform as planned. 
With this in mind, we let c be an additional design variable to be used for decoupling only. 
Optimal decoupling can be pursued by finding c in C(s, c) such that a stochastic norm of 
the off-diagonal terms of C(s, c)G(s, p) is minimized in the frequency range of interest. This 
problem will be considered in Section 4.4. 

3 ROBUST STABILITY 

In this framework, we say that a LTI system with structured uncertainty is robust stable if 
all its poles reside in the left hand side of the complex plane for all possible values of the 
uncertain parameters, i.e the system is asymptotically stable Vp G A. The requirement (5) 
listed above implies the robust stability of the closed loop response. 

Denote with s a vector formed by the poles of a random transfer function. Stability 
conditions can be written as 

P{Re{s(p, k)} < 0} = e (5) 

where P{A} is the probability of occurrence of the event A and e = {e^} satisfies e* G [0, 1] 
for i = 1, . . . , n. Defining q = Re{s}, these conditions can be alternatively written as 

Fq(0)>1-6 (6) 

where F M (m) is the cumulative distribution function (CDF) of m at m = m. The Routh- 
Hurwitz (R-H) test applied to the corresponding characteristic equation leads to the vector 
inequality r = {r*} > 0, where r* is a R-H determinant. In this context, stability conditions 
can be written as 

Fr{ 0) < e (7) 

where e = {e*} satisfies G [0, 1] for i = 1, . . . n. Robust stability is attained iff e = e = 0. 
In this paper robust stability will be enforced via Equation (7). 
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4 CONTROL DESIGN 


Formulations that involve the solution to optimization problems in the frequency domain 
are proposed below. In each problem, a metric of the controlled response, denoted with y, 
is used to build the cost function. Let T(s, p,k) and T(s) be complex variables associated 
to the current and desired designs respectively and 0(s, p, k) = 0 the closed loop stability 
equation. For SISO systems, take T = L, T = L, and 0 = 1 + L. For MIMO systems take 
T = T, f = f and 0 = det{I + L}. 

4.1 Minimizing the Likelihood of Instability 

Let yi be the ith inequality constraint for closed loop stability, i.e. the ith component of 
Equation (6) or (7), applied to the closed loop characteristic equation. Hence, f Yi [yi) is the 
PDF of the corresponding random variable. Designs that maximize robust stability can be 
found by solving the following unconstrained optimization problem: 


k* = argmin 


EfV.M 

i — 1 


(8) 


This problem will be referred to as PI and its cost as J\. The minimization of J\ shapes the 
PDFs of the R-H determinants such that their excursion into the negative part of the r-axes 
is minimized. Robust stability is achieved when m/ k ,{Ji} = 0. Notice that robust stability 
might not be possible for the chosen control structure iF(s,k). 

4.2 Minimizing the Variability About a Target Response 

Assume that T is a deterministic target transfer function whose magnitude y satisfies the 
control requirements (1-4). Let y = \T(s, p, k)|. Hence, f Y (y,u>) is the PDF of the corre- 
sponding random process. Designs that improve the control performance by reducing the 
variability of the response can be pursued by solving the following constrained optimization 
problem 


k* = argmin 


'/ [ 


[y(co) - y(u>, p, k)] 2 f Y (y,u ) g(u) dy duj 


(9) 


subject to the set of stability constraints, i.e. either Equation (6) or (7), for a given e or e, 
prescribed in advance. In this equation, 4>(o;) is the support of y and g(w) is a weighting 
function used to favor particular frequency ranges. This optimization problem will be referred 
to as P2 and its cost as J 2 . 

The minimization of J 2 concentrates the entire family of PDFs of the output about the 
target. Not only the expected value of y is moved towards y but also the variability around 
it is reduced. Notice that this treatment involves shaping the whole PDF, including all its 
moments and corresponding support. While the open loop requirements (1-4) are pursued 
by manipulating the cost directly, the closed loop requirement (5) is enforced via inequality 
constraints. 

The main advantage of this formulation is the need for only the first two order moments 
of y to calculate J 2 . Notice however that this approach is over-restrictive because it penalizes 
solutions that do not necessarily violate the control requirements. 
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4.3 Minimizing the Likelihood of Unfavorable Performances 

Let y = \T(s, p, k)|. Hence, F Y (y,uo) is the CDF of the corresponding random process. 
Define the non-negative functions hi(uo) in uo G [0, cc?i] and h 2 (uo) in uo G [uo 2 , oo) such that 
a design satisfying y < h\(u) violates the requirements (1) and (3) and a design satisfying 
y > h 2 {uo) violates the requirements (2) and (4). Clearly, hi and h 2 are envelopes for the 
undesired low and high frequency behavior respectively. Excursions of y below hi or above h 2 
are unfavorable from the performance point of view. In this context, designs that improve the 
control performance can be found by solving the following constrained optimization problem 


Fy(hi(uo), uo) g(co) duo + [1 — Fy(h 2 (co),co)\ g{uo) duo 


subject to the set of closed loop stability constraints for a a given probability of instability. 
This problem will be referred as P3 and its cost as J 3 . 

The minimization of J3 reduces the probability of violating requirements (1-4) directly. 
In contrast to J 2 , J 3 does not penalize satisfactory solutions. Notice however, that the 
calculation of the some of the first order moments of yi is insufficient to properly estimate 
J 3 . At this point, let’s us define P ex (u) as Fy(hi(uo)) for uo G [0,aq], as 1 — F Y (h 2 (u > )) for 
uo G [uo 2 , 00 ) and as zero otherwise. Hence, J3 can be obtained by integrating P ex {uo)g{uo). 

4.4 Maximizing the Decoupling Effect of the Compensator 

The probabilistic description of p makes the components of the rq x n 2 matrix D = CG = 
{djk}, random variables. Perfect decoupling is achieved if the random processes of the off- 
diagonal terms are Dirac deltas at zero. This is clearly impossible unless G is known exactly. 
Let yjk = \djk(s, p, c)|. Hence, F Y (y,co) is the CDF of the corresponding random process. 
Providing that the decoupler is causal and stable, the compensator can be designed by 
solving: 


c = argmin 


(u,p, c) 2 f Y (y jk , co)g(uj) dy jk duo 


where fl is the frequency range of interest and <&j k {uo) is the support of f Yjk (yj k ,Lo) at uo. 
This problem will be referred as P4 and the corresponding cost as J4. Constraints on the 
robust stability of C can also be imposed. 

5 EXAMPLES 

In this section, the formulations proposed above are evaluated using problems that admit 
closed form expressions for the PDFs of y and r ? > This allows us to have a good understanding 
of the problem at hand by avoiding numerical errors caused by sampling and asymptotic 
approximations. 

Notice that the set of solutions of the PI problem defines the domain where the solutions 
to the P2 and P3 problems exist. Then, a solution to the PI problem can be used as initial 
condition in the numerical search for the P2 and P3 solutions. In all the examples, p = {a} 
and a G A = [ a~,a + ], / 4 (a) are given. 
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5.1 SISO System 

Consider the plant and the PID controller 


G(s,p) 
K(s, k) 


6(1 + ds) 

(l + CS^j^C + CLS^j 

k\ + k2S + k%s^ 
s 


( 12 ) 

(13) 


where k = {ki, k 2 , k 3}. The magnitude of the resulting loop transfer function is given by 


V = \ L \ 



UJ 


h 


y/e 2 + (a ;a ) 2 

f (l + (du>) 2 )[(k 2 - k 3 io 2 ) 2 + (kito) 2 ] 
(cui) 2 + 1 


(14) 


Notice that the effect of the uncertainty is not noticeable at low frequencies. The random 
process of the output is given by 

f Y {y,uj) = [f A (d) + f A (-d )] (15) 

where a = (l/o;)-y/(h/y) 2 — e 2 . The support of y, i.e. y G [y~,y + ], is bounded by 


V (a 1 ) = min{y(a + ),y(a )}, y + (cj) 


h/\e\ if 0 G A 

max{y(a + ), y(a~)} otherwise 


The characteristic equation for closed loop stability is 0(s, p, k) = as 3 + /3s 2 + js + 8 = 0, 
where a = ac + bdk 3 , f3 = a + ce + b(dki + k 3 ), 7 = e + b(ki +dk 2) and 8 = bk 2 . The R-H test 
leads to the inequalities 77 > 0 for * = 1, 2, 3 where rq = (3 /a, r 2 = — aS and r 3 = 8 /a. In 

general, the constraints can be written as r* = (aQ + Vi)/(aQ + k*) > 0 for i = 1, 2, 3. The 
expressions for Q, r/i, Q and which are non-linear functions of the components of k, are 
omitted due to their extension. The corresponding CDFs are given by 


{ 1 — sign{m,j}[F A (aj) — F A (d J )\ if r ; > Cj/£j an d rrij Q 0 

sign{rrij}[F A (dj) — F A (dj)] if rj < Q/Q and rrij Y 0 (16) 

H{rj - ( aQ + rji)/(aQ + «,-)} otherwise 


( F A (a 2 ) if C2 > 0 

F R2 (r 2 ) = l 1 - F A (d 2 ) if C2 < 0 

y H{r 2 — r/z} otherwise 


(17) 


where rrij = Qkj - rjjQ, a 3 = -Kj/Q for j = 1, 3; d t = (niu - rn)/(Q - nQ) for i = 1, 2, 3 
and H{-} is the Heaviside function. Differentiation leads to 


f K|/a (%)/(0 - r£jY if rrij ± 0 

\ 8{r 3 - (aQ + rji)/ (aQ + kQ)} otherwise 


(18) 
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f ( \ _ f Ia(r 2VIC2I if C2 7^ 0 

JR 2 \ 2 ; | S{r 2 — r/ 2 } otherwise 


(19) 


where 5 { • } is the Dirac function. Notice that a combination of the gains can eliminate the 
randomness of the constraint making it purely deterministic. This has also been observed 
in [2]. The following observations are made: (i) when dj G A, some designs lead to /^(r*) 
with unbounded support, (ii) unbounded (rj are obtained where rn : j(k. a) = 0, i.e. 3 r* 
such that lim ri ^ > .fifRi( r i) = 00 . 

The problems PI, P2 and P3 can now be solved numerically using Equations (15), (16) 
and (17). In the results that follow, the values b = 3, c — 10, d = —2, e = 2/5, e = 0 and 
the functions y = |3(1 — 2s)/(20s(2s + l)(s/3 + 1))|, hi(u) = 1 for co G [10 -4 , 10 -2 ] and 
h, 2 (to) = 3/(4u;) for lu G [0.5, 10 2 ] are assumed. The improper integrals used to build the cost 
functions will be approximated in accordance with these intervals, i.e. uu G [10 -4 , 10 2 ]. 

First, assume that a is a beta random variable with parameters {2, 10} and support 
A = [0, 2], In this setup, the uncontrolled system is robust stable. The corresponding bode 
diagram is shown in figure 3. The system’s response at high frequencies is clearly inadequate, 
then some control action is needed. The P2-solution leads to the diagram shown in figure 4. 
Here, g(uo) has been used to favor the performance at low frequencies. The corresponding 
PDFs for the R-H determinants are shown in figure 5. Notice that three supports extend 
towards positive infinity. The same information is shown in figures 6 and 7 for a cost function 
that favors the performance at high frequencies. 

Figure 8 shows the stochastic bode diagram for a robust stable design. Notice the un- 
desired performance regions created by h\(iu) and /^(w) and the corresponding P ex . This 
design is used as initial condition for solving the P3-problem. The corresponding solution 
is shown in figure 9, where e = 0.005(1, 1, 1}. The PDFs of the R-H determinants and the 
stochastic Nyquist diagram are shown in figures 10 and 11. Due to e / 0 s = - lis encircled 
by a portion of the PDF of L and a tail of /# 2 (r 2 ) is in the negative part of the r 2 -axis. In 
all three optimal solutions the control is strictly proper, i.e. = 0. 

Now, assume that a is beta distributed with parameters (2, 2} and support A = [—0.5, 1]. 
The corresponding uncontrolled plant is nominally stable, i.e. G(s, F[p]), but robust unsta- 
ble. This can be seen from figure 12 where the PDFs of the R-H determinants are shown. 
It is interesting to notice that the support of those distributions is given by sets with the 
form S = (— 00 , T] U [$, 00 ), where T and $ are finite values that depend on a~ and a + . 
This occurs due to singularities in r* and A caused by vanishing denominators. This can be 
explained by thinking of a as a moving parameter within A. When a approaches zero from 
the right, a pole moves towards positive infinity and the system is stable . At a = 0, such 
a pole disappears at infinity, i.e. the relative degree of the plant decreases by one, and the 
system stability is now given by a degenerated characteristic equation. When a is decreased 
further, such a pole moves from minus infinity towards zero making the response unstable. 
This analysis explains the non-connected support of the bi-modal PDF. 

The first task is to find k for robust stabilization. The stochastic bode diagram and the 
PDFs for the R-H determinants corresponding to a solution of the PI problem are shown 
in figure 13 and 14. It is observed that: (i) there is a substantial difference between the 
nominal and the expected response at high frequencies in both magnitude and phase, (ii) 



the uncertainty on the sign of a makes the phase plot to spread over a range 180 degrees wide 
at high frequencies, (iii) robust stability cannot be achieved without the derivative action 
and (iv) rq behaves as a deterministic constraint at the optimum, i.e. /^(rq) ~ <5(rq — 10 -3 ). 
Notice that observation (iii) precludes the search for designs on the PID control structure 
that lead to satisfactory performances at high frequency. The P 1-solution leads to figures 
15, 16 and 17. 

5.2 MIMO System 

The following model is derived by studying the angular velocity control of a satellite spinning 
about one of its principal axes [12] 


G(s,p) 


1 

s 2 + a 2 


s — cl 2 a(s T 1) 

— a(s “l - l) S — CL 2 


( 20 ) 


The plant has a pair of poles at s — ±aj so it needs to be stabilized. A diagonal PID control 
structure without compensator is considered, i.e. k = {£q, k 2 , & 3 }, i^(s,k) = (Aq + k 2 /s + 
k 3 s)I, C(s, c) = I. Equation (4) leads to 


v = r = p 


/a/i t ^ 

\a CO I 


+ (1 - p) 


[L\j 1 "cL 2 

la — u)\ 


(21) 


where p = yqc For the sa k e 0 f simplicity in the notation we will 
use a and a to denote the minimum and maximum singular values of L respectively. The 
corresponding random processes are given by 


fa(a,uj) = sign{n 
/s(d,w) = sign{y 


a} 


r / \ da^ \ £ ( \ d^A 


( 22 ) 

(23) 


Denote with a* the real 
y* = y(a*i) and y\ = y(a* 2 ) 
the random process of y is 


( A(a 2 )a 2 - 


fr{y,u) 


< Ia (fil) Oi - 
/ a ( o 4)04 + 
, /a(^ l)hl + 


—a; a 2 =F P\/— 2 (1 + a ’ 2 ) — M 2 
p 2 — a 2 

LOG 2 A a ’ 2 ) — (“ 2 


roots of Equation (21) such that a* < dj for j > i. Denote with 
the extrema of Equation (21) such that y\ < y\. If a = da/dy, 
given by 


f A {d i)ai 

f a (d 2 ) d 2 

/a(o 2)^2 — _/U(d i)hi — fA(d 3 )a 3 
fA(d 3 )a 3 - / 4 (a 2 )a 2 - / A (a 4 )a 4 


if (yl<y < £< y 2 ) or 
(y{ <y< y 2 < aO 
if y{ < £ < y < y 2 
if y* < y 2 < y < n 

otherwise 


(24) 
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Analytical expressions for the as and the y * s were derived but omitted due to their length. 
The support of y G [ y~,y + ] is bounded by 

and a} 0 A for i ^ j 
and G A 


if |u;| G A 
otherwise 

Notice that fy(y,u) is unbounded at y\ and y\ and its support is unbounded from above 
for all possible resonant frequencies, i.e. Mu G A. 

The closed loop stability equation is 0(s, p, k) = as 4 + /3s 3 + qs 2 +x^>s + 0 = 0, where a = 
1 +2&3 + (1 , (5 = 2[ — bk^ -{-k 2 (l + k% FbkF )} , q = k2~\-2ki(^l-\-k^ -\-b(l — 2k2~\~k2^ -\-2k\k^) , 
i[j = 2k\[b(k 2 — 1) + k 2 ] and 6 = (1 + b)kf where b = a 2 . The R-H test leads to the non-trivial 
inequalities r\ = (5 > 0, r 2 = /?q — > 0 and r 3 = /?q?/> — aA 2 — 9/3 2 > 0. These constraints 

can be written as r* = + £j& 2 + (}& + Vi for i = 1, 2, 3. The corresponding CDFs are given 

by 


V Vi 

min{y(a + ),y(a )} 


if a* G A 
if a{ G A 
otherwise 


, s _ j max{y} = 00 
' - 1 max {y(a + ) , y{a~)} 


F Ri( r i) 


' F B (b 1) 

1 -Fsih) 

F B (bs) - F B (b 2 ) - F B (h) 

< 1 + F B (b 2 ) — F B (bi) — F B (b 3 ) 
F B ih)-F B {h) 

1 — F B (b 2 ) + F B (bi) 

, # !>* - ??*} 


if (/c* > 0, n = 1) or (/q = 0, = 0, £i > 0) 

if (fq < 0, n = 1) or (/q = 0, = 0, Ci < 0) 

if /q > 0 and n = 3 

if /q < 0 and n = 3 

if /q = 0 and n = 2 and £ > 0 

if /q = 0 and n = 2 and £ < 0 

otherwise 


(25) 


where i = 1,2, 3, F B (b) = Fa^VF) — Fa{—VF) and n is the number of real roots of Q(b) = 
K ib 3 + iib 2 + Qb + r]i — ri = 0, such that Q(F ) = 0 and \ < bj for j > i. Notice that the 
coefficients of this polynomial are non-linear functions of the components of k. Analytical 
expressions for the bs were derived but omitted here due to their length. Expressions for the 
PDFs of the constraints can be derived by differentiating Equation (25). 

The problems PI, P2 and P3 can now be solved numerically using Equations (24) and 
(25). In the results that follow, a is a beta random variable with parameters {6,6}, A = 
[0.5, 1.5], e = 0 and y, hi(u) and h 2 (u) are as before. Figures 18 and 19 correspond to the 
P2-solution. In this case, a strictly proper control leads to a robust stable system in which r 2 
is purely deterministic. The uncertainty on the natural frequency is clearly reflected in the 
support of fy(y, u). Notice that in spite of this, E[y] is bounded. Figure 20 shows the section 
of the k-domain where k% = 0 for A = [0.5, 10]. The support of the corresponding constraints 
is filled with dotted lines. Notice that the bounds of such supports are not always given by 
the deterministic constraints evaluated at the extreme values of A, i.e. r t {lr) / and 
Vi(b + ) 7^ r+. By constraint envelopes we mean a deterministic set of inequalities satisfying 
the stability conditions for e = e = e = 0 . Having this information at hand will considerably 
reduce the computational demands of searching for the optimum. The admissible domain, 
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where robust solution exists, is colored in gray. Such region is bounded by the tightest set 
of constraint envelopes. Figure 21 corresponds to the P3-solution, whose features are very 
similar to figure 19. 


5.3 Decoupling 

Below, a compensator for the plant given in Equation (20) is designed by solving the P4- 
problem given the configuration shown in figure 2. The compensator will be first designed 
via generalized decoupling. Define G(s, E[p]) as the nominal plant. If the desired uncoupled 
behavior is given by Diag{G(s, E[p])}, the resulting compensator is 

(£[a]-s ) 2 £[a](l + s)(£[a] 2 -s) ' 

-E[a\{ 1 + s)(E[a] 2 - s ) (E[a] - s ) 2 

Notice that the compensator is Bounded-Input-Bounded-Output (BIBO) stable and that 
E[a] acts as a design variable, i.e. c = {c} = {E[a\}. Few manipulations lead to 

„ _ — a 2 (c 2 — s) — ac( 1 + s ) 2 — s(s — c 2 ) — (a — c)(l + s)(ac + s ) 

° (a — c)(l + s)(ac + s ) — a 2 (c 2 — s) — ac( 1 + s ) 2 — s(s — c 2 ) 

(26) 


C{s,c) 


(s 2 + E[a] 2 )(l + E[a] 2 ) 


where D(s, c, p) = CG and a = (c 2 — s)/[(s 2 + c 2 )(l + c 2 )(a 2 -I- s 2 )]. Notice that if a is 
deterministic, i.e. / 4 (a) = 8 (a — E[a\), and c = E[a] perfect decoupling is achieved, i.e. 
D(s, c, p) = Diag{G(s, E[p])}. The presence of uncertainty prevents this to happen. In this 
example the compensator is based on the inverse of the plant, what give us as many design 
variables as uncertain parameters. This, however, does not have to be the case. 

The structure of the matrix (26) allows us to define an equivalent metric for the off- 
diagonal error as follows 


_ g\a-c\ rrr - i — 5 

y = y 12 = 2/21 = 7-^ 5- vrc 2 + lo 2 

\a z — uo z 


(27) 


Denote with a* the real roots of Equation (27) provided that a* < dj for j > i. Let y* be 
one of the n extrema of y, i.e. y* = y(a*) such that dy / da = 0 at a = a|, and y* < y* for 
i<j. If a = da/dy, the random process of y is given by 


fr(y,u) 


/ a (« 2)02 - / a ( o i)di 
f A {di)di - /A(a 2 )a 2 

/a(o 4)^4 + /n(d 2)^2 — fA(d i)ai — f A (a 3)03 
f A {d i)di + f A (a 3 )d 3 - f A (d 2 )d 2 - f A (a 4 )a 4 


if C'i holds 
if C 2 holds 
if C 3 holds 
otherwise 


(28) 


where the C'i, C 2 and C 3 are 

( {n = 1 and [(|c| > w, y\ > y > y) or (|c| >u, y > y* > y) or (|c| < w, y\ > y)]} or 
C'i = < {n = 3 and \{y* 2 > y > y\ > y) or (y£ > y > y > y* 2 ) or (2/3 > y > y 2 > y* > y) or 
[ (y*i >y>y)ox(y>yl>y>yl) or (y > y 3 > y 2 > y{ > y)]} 


{ {n = 1 and (y\> y > y)} or 

{n = 3 and [(y| > y > y 2 > y > y\) or (y% > y > y > y* 2 ) or 

(2/3 > y > vl > y{ > v) or (y* > y > y)]} 
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{ {n = 1 and [(|c| >u, y > y\ > y) or (|c| > oj, y > y > y{) or (|c| <u, y> y)]} or 
{n = 3 and [(y > y* z > y% > y > y\) or (y* z > y > y^ > y > y > y\) or (y > y > y* z ) 
or (y > yt > y > y > y$) or (y > y* z > y\ > y\ > y) or (y z >y 2 >y> y[ > y)]} 

and y = \c\y. Analytical expressions for aj, a* and y* can be found. The support of y is 
bounded by 


if c e A 

if a] e A and c ^ A 
if al € A and{a, a^} ^ A and n = 3 
otherwise 

if |w| € A 

if a* 2 E A and \uj\ £ A and n = 3 
^ max{y(a + ), y(a~)} otherwise 

Problem P4 can now be solved numerically using Equation (28). In the results that 
follow the same data used in Example 5.2 is assumed and Q 6 [0, oo). Figure 22 shows the 
magnitude of the off-diagonal terms of D versus frequency for c = E[a] = 1. This is the 
resulting behavior of a conventional decoupling design, where perfect decoupling is achieved 
if the plant is known exactly. Notice that y~ and y(E[a\) coincide. This figure indicates 
significant discrepancies between deterministic and stochastic results. 

The P4-solution leads to c* ^ E[a\. Figure 23 shows a slight improvement in the decou- 
pling at low-frequencies. Compensators with more design variables and/or different struc- 
tures can achieve better results. Note that: (i) the compensator’s structure does not have 
to be based on the inverse of G, (ii) the decoupling of non-square plants can be studied by 
the same means. 

6 DISCUSSION 

In the cases where the uncertain parameters represent a set of changing operating conditions 
this methodology leads to valuable solutions that otherwise are very difficult to find. In 
problems where the parameters are physical quantities modeled as uncertain variables as a 
result of our ignorance on their real value, the solutions can be improved further by refining 
and even eliminating their probabilistic description using additional information. Adaptive 
control and model predictive techniques can be used in this respect. 

The ignorance on the actual implemented values of the solution can be studied in advance 
by modeling physical design quantities as additional random variables. We refer to this 
type of uncertainty as input uncertainty . This can be achieved by using variables that 
parameterize prescribed PDFs as design variables in the optimization problem. Notice that 
while the physical quantity is deterministic, its modeling is stochastic and its manipulation 
is done via deterministic variables. For example, if a gain is modeled as a Gaussian random 
variable with a fixed relation between its first two order moments -say, due to limitations 
in the precision of the solution to be implemented- optimal solutions can be found by using 
a single deterministic variable in the optimization problem. Designs that carry out this 
practice will lead to satisfactory performances for the entire range of gain values within the 


{ min{y} = 0 

i 

mm{y(a + ),y(a )} 
( max{v\ = oo 

?/» = \ y*2 
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support of the corresponding PDFs. Notice that in spite of using the same modeling for 
both parameter and input uncertainty, their physical significance is very different. 

The behavior of the probabilistic constraints and their envelopes in the k-space must be 
explored. Such envelopes cannot, in general, be easily built (providing they exist). This 
fact prevents us from relaxing the computational demands of the optimization problem. The 
tightest constraint envelopes, given by the set that maximizes the size of the robust stable 
region in the k-space, can be obtained by propagating A at a high computational expense. 
On the other hand, over-constraining envelopes might be easier to build at the expense of 
introducing conservatism into the problem. This is especially valid in the cases where at least 
one of the constraints is active (in the probabilistic sense) at the optimum. Circumventing 
the probabilistic nature of the constraints leads to the loss of the sensitivity of the probability 
of instability to changes in k. This information is very valuable from the design point of view. 
This feature posts a formidable trade-off between numerical convenience and mathematical 
significance. 

The control design of systems with multiple uncertain parameters lead to additional 
challenges. In such problems, the calculation of moments and probabilities have to be done, 
almost exclusively, by means of sampling techniques and asymptotic approximations. While 
sampling based techniques are suitable for calculating first order moments and large probabil- 
ities; asymptotic approximations are reliable for small probabilities. Because the probability 
of occurrence of the metrics used in the optimization problem can fluctuate considerably, 
algorithms that dynamically discriminate between these tools are desirable. In addition, the 
repeated evaluation of the cost function suggest the need for adjusting the trade off between 
accuracy and computational cost dynamically. These issues are currently being investigated. 

7 CONCLUSIONS 

In this paper, a stochastic approach to the control synthesis of systems with structured 
uncertainty is introduced. This strategy allows us to quantify the probability of instability 
and/or unfavorable performance such that robust optimality can be pursued. Formulations 
for the stochastic optimization of SISO and MIMO control systems are proposed taking 
into account decoupling, performance and stability considerations. Examples that admit an 
exhaustive analysis of the formulation and its solution are used to validate the approach. 
The application of this methodology to large scale problems, where the use of reliable and 
inexpensive numerical tools is crucial, is currently being addressed. 
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Figure 1: Block diagram of one degree-of- freedom SISO control system 



Figure 2: Block diagram of one degree- of- freedom MIMO control system 



icr 2 io"' io° io' io z 


Figure 3: Stochastic Bode diagram for the SISO uncontrolled system. The PDFs of the 
magnitude and phase are non-zero within the solid lines shown. The nominal plant |G(i£[p])| 
and its expected value F[|G ! (p)|] are indicated with dash-dotted and dashed lines respectively. 
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Figure 4: Stochastic Bode diagram of a design that favors the performance at low frequencies. 
The PDFs are non-zero within the solid lines shown , i.e. y~ < y < y + . The nominal 
plant y(E[j>]) and its expected value E[y( p)] are indicated with dash-dotted and dashed lines 
respectively. The curve of the target function y is shown using a line- circle pattern. The 
same line conventions are used in the bottom for <j>(L). 



Figure 5: PDFs of the R-H determinants for the optimal solution, f^iri), ff. 2 (r^) and 
f Rz (r 3 ) are show using dashed, dash-dotted and solid lines respectively. 
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Figure 8: Top: upper view of fy(y, of) for a robust stable design. Conventions used previously 
apply. Dotted lines are used to show the boundary of the sets defined by y < hi(u) and 
V > h 2 (to). Bottom: corresponding P ex . 
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Figure 9: Top: upper view of fy(y,u}) corresponding to the optimal solution. Conventions 
used previously apply. 
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Figure 10: PDFs of the R-H determinants for the optimal solution. Conventions used previ- 
ously apply. 



Figure 11: Stochastic Nyquist plot corresponding to the optimal solution. Crosses centered at 
the expected value are shown for discrete frequency values. Dotted lines are used to indicate 
the support of y in the radial and tangential directions. A dash-dotted line and a dashed line 
are used to show the nominal and expected Nyquist plots respectively. 
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y(E[p]), E[y(p)] 



Optimal solution to the P2 problem. y(E[p]), E[y( p)], y~ and y + are indicated 
tted, dashed and solid lines respectively. The curve of the target function y is 
a line-circle pattern. 





Figure 20: Section of the k-domain for k 3 = 0. The support of the constraints is filled with 
dotted lines and the functions rj(a + ) and rj(a _ ) are shown using solid lines. Robust stable 
regions, where the optimal resides, are colored in gray. 
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Figure 21: Top: upper view of fy(y) for a robust stable design. Conventions used previously 
apply. Bottom: corresponding P ex . 
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Figure 22: Results for c = E[a]. Top: magnitude of the off-diagonal term of B. The PDFs 
are non-zero within the solid lines, i.e. y~ < y < y + . The nominal and expected behaviors, 
namely y(E[ p]), E[y( p)] are marked with dash-dotted and dashed lines respectively. In this 
case, y~ and y(E[ p]) coincide. Bottom: probability of the norm of the off-diagonal term to 
exceed yu m = 0.05. 



Figure 23: Results to the Pf-problem. Conventions and values used in the previous figure 
apply. 
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